Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
mkdir ~/EVALUATION
cd ~/EVALUATION
mkdir RAW_DATA REFERENCES QC CLEANING MAPPING
module load sra-tools/2.10.3
fasterq-dump --version
# "fasterq-dump" version 2.10.3
srun --cpus-per-task 8 fasterq-dump -S -p SRR10390685 --outdir RAW_DATA --threads 8
cd RAW_DATA
gzip SRR10390685_1.fastq
gzip SRR10390685_2.fastq
cd -
cd ~/EVALUATION/RAW_DATA
module load seqkit/0.14.0
seqkit stat *.fastq.gz > seqkit.tsv
head seqkit.tsv
# file format type num_seqs sum_len min_len avg_len max_len
# SRR10390685_1.fastq.gz FASTQ DNA 7,066,055 1,056,334,498 35 149.5 151
# SRR10390685_2.fastq.gz FASTQ DNA 7,066,055 1,062,807,718 130 150.4 151
# other solution but less generic
# zcat SRR10390685_1.fastq.gz | echo $((`wc -l`/4))
# zcat SRR10390685_2.fastq.gz | echo $((`wc -l`/4))
Les fichiers FASTQ contiennent 7 066 055 reads.
cd ~/EVALUATION/REFERENCES
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
gzip -d GCF_000009045.1_ASM904v1_genomic.fna.gz
cd ~/EVALUATION/REFERENCES
module load seqkit/0.14.0
seqkit version
# seqkit v0.14.0
seqkit stat GCF_000009045.1_ASM904v1_genomic.fna
# file format type num_seqs sum_len min_len avg_len max_len
# GCF_000009045.1_ASM904v1_genomic.fna FASTA DNA 1 4,215,606 4,215,606 4,215,606 4,215,606
La taille de ce génome est de 4 215 606 paires de bases.
cd ~/EVALUATION/REFERENCES
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
zgrep -v "^#" GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '($3 == "gene")' |wc -l
# 4448
4 448 gènes sont recensés dans le fichier d’annotation.
cd ~/EVALUATION/QC
module load fastqc/0.11.9
fastqc --version
# FastQC v0.11.9
for i in ../RAW_DATA/*.fastq.gz ; do srun --cpus-per-task 8 fastqc $i -o . -t 8 ; done
module load multiqc/1.9
multiqc --version
# multiqc, version 1.9
srun multiqc -d . -o .
La qualité des bases me paraît satisfaisante car la qualité moyenne des bases est supérieure à 30 sur toute la longgueur des reads comme le montre le graphique Sequence Quality Histograms ou le graphique Per Sequence Quality Scores. De plus, il y a très peu de bases indéterminées (N).
Lien vers le rapport MulitQC
Oui, car la distribution de taille des reads (Sequence Length Distribution) montre que certains reads sont plus petits que 151, la longueur attendue. Néanmoins, cela représente très peu de reads donc ce n’est pas vraiment problématique. Il serait intéressant de trouver l’information des filtres appliqués (via la publication pour des données publiques ou via la plateforme de séquençage pour des données non publiées.
Formule à appliquer :
(length all R1 + length all R2) / Taille du génome
soit ici:
(1056334498+1062807718)/4215606
La profondeur de séquençage est de : 502.6898187 X.
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp qui vous semblent adéquats et justifiez-les.
cd ~/EVALUATION/CLEANING
module load fastp/0.20.0
fastp --version
# fastp 0.20.0
srun --cpus-per-task 8 fastp --in1 ../RAW_DATA/SRR10390685_1.fastq.gz --in2 ../RAW_DATA/SRR10390685_2.fastq.gz --out1 SRR10390685_1.fastq.gz --out2 SRR10390685_2.fastq.gz --html fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json fastp.json
seqkit stat SRR10390685_[12].fastq.gz > seqkit.txt
head seqkit.txt
# file format type num_seqs sum_len min_len avg_len max_len
# SRR10390685_1.fastq.gz FASTQ DNA 6,777,048 996,891,051 100 147.1 151
# SRR10390685_2.fastq.gz FASTQ DNA 6,777,048 990,442,597 100 146.1 151
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| –cut_mean_quality | 30 | pour un score moyen dans la fenêtre glissante > 30 |
| –cut_window_size | 8 | pour une taille de fenêtre glissante de 8 |
| –length_required | 100 | pour ne garder que les reads de taille > 100 |
| –cut_tail | pour faire partir la fenêtre de l’extrémité 3’ du read |
Ces paramètres ont permis de conserver 6 777 048 reads pairés, soit une perte de 4.0900757 % des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [3] et samtools [4].
cd ~/EVALUATION/MAPPING
module load samtools/1.10
samtools --version
# samtools 1.10
# Using htslib 1.10.2
module load bwa/0.7.17
bwa
# Version: 0.7.17-r1188
srun bwa index ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.fna
srun --cpus-per-task=4 bwa mem ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.fna ../CLEANING/SRR10390685_1.fastq.gz ../CLEANING/SRR10390685_2.fastq.gz -t 3 | samtools view -hbS - > SRR10390685.bam
srun samtools flagstat SRR10390685.bam > SRR10390685.bam.flagstat
srun samtools sort SRR10390685.bam -o SRR10390685_sorted.bam
srun samtools index SRR10390685_sorted.bam
samtools view -f 4 -c SRR10390685.bam
# 744540
744 540 reads ne sont pas mappés.
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [5]:
Pour répondre à cette question, je choisis de récupérer tous les reads qui chevauchent le gène en faisant l’intersection des reads du fichier BAM avec le fichier GFF3. Ensuite je compte le nombre de reads du fichier BAM généré.
cd ~/EVALUATION/MAPPING
module load bedtools/2.29.2
bedtools --version
# bedtools v2.29.2
zgrep trmNF ../REFERENCES/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > trmNF.gff3
bedtools intersect -a SRR10390685_sorted.bam -b trmNF.gff3 -f 0.5 > SRR10390685_on_trmNF.bam
samtools view -c SRR10390685_on_trmNF.bam
# 2801
2 801 reads chevauchent le gène d’intérêt.
Utilisez IGV [6] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
Pour répondre à cette question, il faut récupérer en local le fichier BAM et son index (.bam.bai) ainsi que le génome de référence, et éventuellement le GFF3 du gène d’intérêt pour observer les parties flanquantes.
IGV : couverture du gène trmNF
tree ~/EVALUATION
EVALUATION
├── CLEANING
│ ├── fastp.html
│ ├── fastp.json
│ ├── seqkit.txt
│ ├── SRR10390685_1.fastq.gz
│ └── SRR10390685_2.fastq.gz
├── MAPPING
│ ├── SRR10390685.bam
│ ├── SRR10390685.bam.flagstat
│ ├── SRR10390685_on_trmNF.bam
│ ├── SRR10390685_sorted.bam
│ ├── SRR10390685_sorted.bam.bai
│ └── trmNF.gff3
├── QC
│ ├── multiqc_data
│ │ ├── multiqc_data.json
│ │ ├── multiqc_fastqc.txt
│ │ ├── multiqc_general_stats.txt
│ │ ├── multiqc.log
│ │ └── multiqc_sources.txt
│ ├── multiqc_report.html
│ ├── SRR10390685_1_fastqc.html
│ ├── SRR10390685_1_fastqc.zip
│ ├── SRR10390685_2_fastqc.html
│ └── SRR10390685_2_fastqc.zip
├── RAW_DATA
│ ├── seqkit.tsv
│ ├── SRR10390685_1.fastq.gz
│ └── SRR10390685_2.fastq.gz
└── REFERENCES
├── GCF_000009045.1_ASM904v1_genomic.fna
├── GCF_000009045.1_ASM904v1_genomic.fna.amb
├── GCF_000009045.1_ASM904v1_genomic.fna.ann
├── GCF_000009045.1_ASM904v1_genomic.fna.bwt
├── GCF_000009045.1_ASM904v1_genomic.fna.pac
├── GCF_000009045.1_ASM904v1_genomic.fna.sa
└── GCF_000009045.1_ASM904v1_genomic.gff.gz
6 directories, 31 files
1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.
4. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.
5. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
6. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France